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We consider the non Abelian sandpile model introduced by Y.-C. Zhang on a two-dimensional 
square lattice. The static and dynamical properties of the model are investigated and compared to 
the Abelian sandpile model of Bak, Tang and Wiesenfeld. A detailed analysis which takes the finite 
size effects into account yields that the exponents of the avalanche probability distribution are the 
same as in the Abelian model. 
PACS number: 05.40. +j 
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I. INTRODUCTION 

The idea that an externally driven physical system 
with many degrees of freedom can be self-organized crit- 
ical (SOC) was introduced few years ago by Bak, Tang 
and Wiesenfeld (BTW) and realized theoretically using 
a stochastic cellular automaton The original BTW 
model belongs to the Abelian sandpile models Here 
the sequence of relaxation processes is described by op- 
erators which satisfy a commutative algebra. This prop- 
erty allows the analytical calculation of some features of 
the system in the steady state A continuous ver- 

sion of this model was introduced by Zhang to study the 
propagation of activated energies Q . In contrast to the 
BTW model the Zhang model is a non Abelian model, 
i.e., the steady state configurations depend on the se- 
quence in which unstable sites are toppled (see and 
references therein). Despite the different microscopic dy- 
namics both models are expected to belong to the same 
universality class (see for instance ||^). Up to now no- 
body has proved this assumption by direct measurements 
of the avalanche exponents on large lattice sizes which 
reduces the finite size effects sufficiently. We consider 
the Zhang model on lattice sizes which are significantly 
larger than those sizes used in previous investigations 
|,[0|-[l2). The energy distribution p{E) which charac- 
terizes the static properties of the model is concentrated 
around z distinct peaks where z is the number of nearest 
neighbors. We show that the peaks are located at mul- 
tiples of and that the height of the peaks grow with 
increasing system size. Numerical simulations of the two 
dimensional square and honeycomb lattices confirm this 
result. We also investigated the avalanche distributions 
on lattice sizes up to i = 2048. A finite size analysis of 
the exponents of the avalanche distributions yields values 
which corresponds to that of the BTW model. 



II. MODEL 



We consider a two dimensional square lattice of linear 
size L. A continuous value Eij > representing the en- 
ergy is associated to each lattice site (i, j). The boundary 
sites are fixed to zero {Eiboundary) = 0) for all times. A 



configuration {Ei j} is stable if Ei j < E^ for all lattice 
sites («,j). For the sake of simplicity we choose in all 
simulations Ec = 1- A quantum of energy S is added to 
a randomly chosen lattice site (i, j), i.e.. 
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In the case that due to this perturbation a site exceeds 
the critical value Ec an activation event will occur and 
the critical site relaxes to zero and the energy is added 
to the next neighbors, i.e.. 
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where z denotes the number of next neighbors. In that 
way the transferred energy may activate the neighboring 
sites and thus an avalanche of relaxation events may take 
place. Energy may leave the system only at the bound- 
ary. 

In our simulations we use various values of the input 
energies out of the interval 5 g]0, Ec] . In the case of 5 — > 
all lattice sites grow parallel. In order to implement this 
different perturbation process one has to find the site 
with the largest energy Emax and then increment all sites 
by Ec — E,nax ■ In this case the Zhang model is identical 
with the conservative limit of the "spring block" model 
of Christensen and Olami |Q . 

The concept of self-organized critical systems refers to 
driven systems which organize themselves into a steady 
state. We consider the average energy 



1 
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to check if the system has reached the steady state. Start- 
ing with an empty lattice we consider the growth of 
the pile. In the beginning all sites are subcritical, i.e., 
Ei_j <^ Ec and no toppling event occurs. Here, no relax- 
ation process takes place (non avalanche regime) and the 
average energy increases linear in time (see Fig. nf) . With 
further perturbations the average energy is still growing 
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FIG. 1. The average energy {E{t)) as a function of the FIG. 2. The probability distribution p(i5) for different sys- 

rescaled time r = 5L~^t for L < 512 and various values of 5. tem sizes. The inset displays the scaling plot of the third 

maximum of p{E). 



until one site reaches the critical value Ec- Now the be- 
havior of the system changes and toppling processes oc- 
cur (avalanche regime). After a certain time the average 
energy reaches a constant value (E) which characterizes 
the steady state. 

In Fig. the average energy is plotted as a function 
of the rescaled time t = 6 t. One can see a data 
collapse of all curves corresponding to different values 
of L and 6. Deviations from the collapse occur only at 
the point t « 0.63 where the behavior changes from the 
non avalanche regime to the avalanche regime, character- 
ized by a constant average energy. The avalanche regime 
occurs when the fluctuations of the energies are greater 
than the difference of the critical energy Ec from the av- 
erage energy {E{t)), i.e., when 

ViEHt)} - {E{t))^ > E,~ (Eit)). (5) 

Decreasing S reduces the fluctuations and the critical 
time tends to Tc — 1. Larger system sizes result in a 
decreasing critical time. 

We consider the system for lattice sizes L G 
{64, 128, 256, 512, 1024, 2048} (in the case of 5 = the 
maximum lattice size is L = 1024). Starting with an 
empty lattice the system will be equilibrated after L^S^^ 
perturbations. In order to provide a sufficient statistics 
all measurements are averaged over at least 10^ non zero 
avalanches. 



III. ENERGY DISTRIBUTION 

We measured the energy distribution p{E) in the 
steady state for S G {0, 128^^, 8"-^, 1} and various sys- 
tem sizes L. In Fig. ||the distribution p{E) is plotted for 



different system sizes. The distribution is concentrated 
around four distinct peaks. It was assumed in previ- 
ous works ppl] that the finite spreads of the peaks are 
caused by "intrinsic dynamical fiuctuations" . As one can 
see from Fig. || the peaks grow and the spreads of the 
peaks decrease with increasing system size L. 

We found that the maximum Pmax{E) of each peak 
scales with the system size as 

PmaAE) i^ (6) 

with y « 0.6. Since the distribution p{E) is normal- 
ized we assume that the peaks scale in the horizontal 
direction as L^^ . The location of the maxima of the dis- 
tribution p{E) depend slightly on the system size L. In 
order to produce a scaling plot this drift has to be taken 
into account. In the inset of Fig. ^ we plot L~yp{E) 
as a function of Ly{E — Emax{L)) and get a satisfying 
data collapse. The peaks of the energy distribution p{E) 
grow to infinity and the spread of each peak vanishes for 
L oo. In the case of an infinite system the energy 
distribution p{E) in the steady state is given by 

3 

p(i?) = ^ /,5(i?-i?0, (7) 

where fi denotes the statistical weight and Ei denotes 
the location of the (5-peaks. 

One can calculate the discrete values of the energies Ei 
in the following way [Q: Suppose that the energies are 
already discretized with the allowed values 

Ee{0,Eo, 2Eo, 3Eo, ... , nEo, ...}. (8) 

Then a maximum value of n exists with 
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FIG. 3. The statistical weights fi as a function of the in- 
verse system size for d = 8~^. The values /i(oo) are 
obtained by an extrapolation to the vertical axis. 



FIG. 4. The probability distribution Ps (s) for different sys- 
tem sizes for S — 128~^. The curves for L < 2048 are shifted 
in the downward direction. 



Eq < Ec < (n„ 

The critical energy E — {rimax 
should be equal to i5o, i-e., 



+ l)Eo. (9) 
'1)Eq relaxes and — 



l)Ea 



Eq. 



(10) 



In this way the number of peaks equals the lattice coor- 
dination number rimax + 1 = z. Based on his numerical 
investigations of different lattice types Dfaz-Guilera has 
already proposed this relation | |T^ ]. 

Starting with a stable configuration one perturbs the 
system until one site becomes critical, i.e., one adds 
AE = Ec — Umax Eo on each lattice site (this is correct for 
(5^0). The energy of a given site is now E = uEq + AE. 
The critical site relaxes and is added to the z next 
neighbors of this site. Arguing that the new energy is 
the next allowed energy value E = {n + 1)Eq one gets 
the relation 
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Ec 



z + 1 



(11) 



Note that the discretization of the energies is independent 
of the dimension of the system. The relevant term is the 
lattice coordination number z. This is in contrast to the 
conclusions drawn from previous investigations. These 

TABLE I. Statistical weights of the energy distribution. 
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based only on simulations of square lattices in different 
dimensions d where the coordination number is given by 
z = 2d [|0|. 

We compare Eq. (jl^) with the results obtained from 
computer simulations. In = 2 we found that Eq ~ 
0.3149 for S = 128-\ Eq ~ 0.3153 for 5 = 0, Eq 
0.3145 for S = 8-\ and Eq ~ 0.3140 for (5 = 1 which 
are in good agreement with Eq. (^ij). We also measured 
the energy distribution of a honeycomb lattice in two 
dimensions (z = 3) and found for 5 = 128~^ the value 
Eq fa 0.443 which corresponds very good to the exact 
value Eq — 0.4. Pietronero et al. have investigated the 
rf = 3 Zhang model on a square lattice and found 6 peaks 
in the energy distribution [|ll| . We measured the average 
distance between two peaks from the Fig. 2 of |11| and 
determined in this way Eq = 0.190 which agrees with 
Eq = 0.194 obtained from Eq. (|lj). 

Furthermore, we determined the statistical weights fi 
of the energy distribution [Eq. (0)] . We divided the in- 
terval [0, Ec] in four parts and measured in each part the 
area fi (L) under the curve p{E) for various system sizes 
L. The statistical weights fi are given by an extrapola- 
tion to L — > cx) (see Fig. ^) and the obtained values are 
listed in Table ^. Analogous to the locations of the peaks 
the statistical weights do not depend on the input energy 
S. On the other hand one can see that the values differ 
from those of the BTW model which are known exactly 

Pietronero et al. [g^ introduced a renormalization 
group approach for sandpile models where the density of 
the critical sites determines the fixed point of the renor- 
malization transformation. Here, the density of the crit- 
ical sites corresponds to the statistical weight /s in the 
sense that any perturbation of a coarse grained particle 
(Eq) leads to a relaxation event. Following our results 
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FIG. 5. The values of the exponent Ts as a function of the 
input energy 5 for a fixed system size L. Note that the values 
of the exponent are independent of 5 in the limit 5 <^ Ec. 



FIG. 6. The system size dependence of the exponent Ts for 
(5 = and S = 128~^. The inset displays the determination 
of Too according to Eq. ( p^ (dashed line). 



both models are characterized by different fixed points 
and thus one might expect that both models belong to 
different universality classes. But one has to emphasize 
that this renormalization group approach and its im- 
provement by Ivashkevich js) neglects fluctuations at the 
steady state. Due to this "mean-field-type approxima- 
tion" we think that the different critical densities of 
the Zhang and the BTW model cannot lead to an answer 
of the universality question. 



IV. AVALANCHE DISTRIBUTIONS 

In this section we examine the probability distribution 
of an avalanche of size s, area duration t, and radius 
r, where s denotes the total number of toppled sites, 
is the number of distinct sites which corresponds to the 
area of an avalanche. The duration t of an avalanche is 
equal to the number of update sweeps needed until all 
sites are stable again. The linear size of an avalanche r 
is measured via the radius of gyration of the avalanche 
cluster. In the critical steady state the corresponding 
probability distributions should obey power-law behavior 
characterized by exponents r^, t^, t^, and according 
to 



Ptit) ~ 

Pr{r) - r"^"-. 



(12) 
(13) 
(14) 
(15) 



The distribution Ps{s) is plotted in Fig. | for 5 = 128^^ 
and various system sizes L. All curves fit in the middle 
region to a straight line and the corresponding exponents 
are determined via regression of this region. First we 
investigate whether the exponents depend on the input 
energy 5 and second we examine how the finite system 
size affects the results. Figure |^ shows the exponent 
for L = 256 and for various values of 5. In the limit 
5 Ec — 1 the exponents are independent of the input 
energy. This behavior changes abrupt for (5 > 32 ~^ where 
the exponent displays a complex 5 dependence. In the 
following we focus our attention on the limit S <ti Ec- 

The exponents Ts corresponding to different values L 
and 5 are plotted in Fig. ^. Significant differences be- 
tween the values of the exponents Ts{L,6 — 0) and 
Ts{L,S = 128^^) are caused by the system size only 
and not by the input energy. Both exponents tend to 
Ts « 1.28 with increasing L. In order to determine the 
exact value of the exponent Tj we assume that its system 
size dependence is given by 



TsiL) 



const 



(16) 



We tried several values of x and got the best results for 
X — 1, i.e., the finite size effects are of the relative magni- 
tude of the boundary (~ L~^). In the inset of Fig. ^ the 
exponents Ts{L) are plotted as a function of the inverse 
system size. The exponent Ts is given by an extrapolation 
to L ^ cx) which yields = 1.282 ± 0.01. 

The exponents of the avalanche probability distribu- 
tion of the area and radius are characterized by the same 
finite-size corrections. The exponents corresponding to 
different system sizes are plotted in Fig. 0. Except of the 
deviation for L = 64 in the case of the exponent both 
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FIG. 7. The values of the exponent Td and Tr as a function 
of the inverse system size for S — 128~^. The dashed 
lines correspond to the extrapolation according to Eq. ( p^ ) 



FIG. 8. The probability distribution Pt{t) for a fixed sys- 
tem size. The dashed line corresponds to a power-law with 
the exponent of the BTW model rt = | (see Q). 



exponents depend on the inverse system size (correspond- 
ing to Eq. (p^)). From the extrapolation to the infinite 
system size we obtain the values Td = 1.338 ±0.015 and 
Tr = 1.682 ± 0.018, respectively. The finite size depen- 
dence explains why lower values of the exponents were 
reported in previous works based on numerical simula- 
tions of one system size only (see for instance [0). 

The finite size analysis described above fails in the case 
of the duration exponent rj. Here, the probability dis- 
tribution exhibits a finite curvature which makes it im- 
possible to determine the exponent via regression (see 
Fig. H) . Using a momentum-space analysis of the corre- 
sponding Langevin equations Dfaz-Guilera showed that 
the dynamical exponent of the BTW and Zhang's model 
is given by z = {d + 2)/3 This result allows to deter- 
mine the exponent Tt because the exponents z,Tt, and 
have to fulfill the scaling relation (see for instance ||l7|l ) 



(17) 



1 



Using the above value of Tr and z — ^ for the two- 
dimensional model we obtain the value = 1.512±0.014. 

Recently, it has been shown numerically that the expo- 
nents of the BTW model are consistent with the values 
Tt = I , = I , and Tr — ^ [|l^ . Because of the lack of a 
scaling relation the exact value of Ts is still unknown but 
the authors estimate the value — 1.293 ±0.009. These 
values are in agreement with our results, strongly sug- 
gesting that both models are characterized by the same 
exponents. 

Note that we determined the exponents of the Zhang 
model for the limit S <ti Ec only. The measurements for 
a fixed system size and larger values of the input energy 
5 yield different values of the exponents (see Fig. ||). But 



this does not mean that the exponents of the infinite 
system size depend on S. It is also possible that the 
finite size behavior [Eq. ([l^)] changes outside the limit 
5 <^ Ec- Further work has to be done to examine how 
the finite system size affects the values of the avalanche 
exponents for 5 k. Ec- 



V. CONCLUSIONS 

We have studied numerically the static and dynami- 
cal properties of the non Abelian Zhang model on large 
system sizes. The steady state energy distribution is con- 
centrated around z distinct peaks which are located at 
multiples of where z denotes the lattice coordina- 
tion number. The statistical weights of the peaks are 
independent of the input energy 5 but differ from those 
of the BTW model. A finite size analysis of the avalanche 
probability distributions in the limit 5 <^ Ec yields ex- 
ponents which are in agreement with the values of the 
exponents of the BTW model. Both models belongs to 
the same universality class, i.e., both models displays the 
same large scale behavior, characterized by the avalanche 
exponents. 
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